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The state of art in studying thermodynamic properties of hot and dense nuclear 
matter is reviewed with the special emphasis on the confinement-deconfinement 
transition between hadron matter and quark-gluon plasma. The most popular mod- 
els used for describing deconfinement are analysed, including statistical bootstrap 
models, pure phase models, the model of clustered quarks, and the string-flip po- 
tential model. Predictions of these models are compared with the lattice numerical 
simulations. It is concluded that precursor fluctuation effects must be taken into 
account in order to get a realistic description of deconfinement transition. The ex- 
istence of precursor fluctuations is in line with the dynamical confinement scenario 
and suggests that deconfinement cannot be considered as a transition between pure 
hadron and quark-gluon phases. All this supports the concept of cluster coexistence 
advocated by the authors of this review: Quark-gluon plasma and hadron clusters 
are different quantum states of the same system, so that any statistical model pre- 
tending to treat nuclear matter under extreme conditions must incorporate into itself 
the probability of these different channels. The ways of constructing statistical mod- 
els with plasma-cluster coexistence are discussed and thermodynamic properties of 
such models are analysed. 



1 Introduction 

One of the most intriguing problems of high energy physics is the possibility of 
transforming the nuclear matter composed of hadrons into the phase consisting 
of their fundamental constituents, quarks and gluons. This phase, because of the 
apparent analogy with the electron-ion plasma, is called the quark-gluon plasma. 
The transformation of the hadron matter into the quark-gluon plasma is named 
deconfinement; the inverse process, respectively, being called confinement. The de- 
confinement transition is somewhat similar to the ionization of atoms. The literature 
devoted to this phenomenon is so numerous that, not to overload the list of refer- 
ences, we shall cite here mainly review papers, when these are available, and in 
which the reader can find thousands of references to original works. The specific 
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feature of the present review is that we concentrate on the statistical models of 
strongly interacting systems under extreme conditions, when non-hadronic degrees 
of freedom become important and deconfinement occurs. 

The possibility that quark degrees of freedom can come into play in the pro- 
cess of relativistic nuclear collisions at already achieved accelerator energies was 
advances by A.M.Baldin [1] who predicted and explained the commulative effect as 
a manifestation of the formation in the colliding nuclei of multiquark droplets. 

In order that quark-gluon degrees of freedom would be essential, special con- 
ditions are necessary, like high temperature or density. From quantum chromody- 
namics it is known that at asymptotically high temperature quarks and gluons are 
really deconfined forming the quark-gluon plasma [2-6]. We also know that at zero 
temperature and at the normal density of nuclear matter there is complete confine- 
ment so that only hadrons exist. But where is the intermediate region in which 
quark-gluon degrees of freedom become relevant? The characteristic parameters for 
this region can be estimated as follows. 

Nuclear matter in the normal state has the baryon density Hqb — 0.167/ fm 3 . 
The corresponding normal quark density is po = 3no.B = 0.5/ fm 3 . The character- 
istic quark interaction energy in the normal state can be presented as E = h/r , 
with the interaction time r = a /c , where ao = p ^ i s the mean interquark 
distance in the normal state. Accepting the system of units in which h — c — 1 , 
with the conversion constant he = 197.327 MeV fm , we have 

E = po /3 = 157 MeV. 

Note an interesting relation between this interaction energy and the characteristic 
baryon energy density Eb = ^jv^ob , in which m^ = \{m p + m n ) = 939 MeV is 
the average nucleon mass. Since Eb = 157 MeV/ fm 3 , therefore E = eb fm 3 . 
Thus, the normal nuclear matter should start to decompose being heated up to 
the temperature equal to the quark interaction energy E , that is up to 6 C ~ 
160 MeV . 

If one wishes to destroy the hadron matter by compression, one has to reach a 
density of about the density of quarks inside a nucleon. This characteristic density 
is p c = 3/vjy , where vjy = (47r/3)r^ is the nucleon volume. Taking for the nucleon 
radius r N = 0.9 fm , we get vn = 3 fm 3 . From here 

Pc = 2p = 1/fm 3 . 

This tells us that the hadron matter should start desintegrating being compressed 
up to the density p c ~ 2p • 

As we see, the expected critical temperature and density are fairly low. Such 
conditions certainly existed in the early Universe about 10 -5 -h 10~ 4 sec after the 
Big Bang and are likely to exist in the interior of neutron stars [7]. 

More important is the common belief that such conditions can be created in 
the process of relativistic heavy ion collisions, even with existing accelerators, thus 
opening the path for experimental observation of the quark-gluon plasma in the 
laboratory [8-14]. Analogous relativistic nuclear collisions can also be studied in 
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cosmic ray experiments [8]. The density of matter inside a fireball formed by two 
collided ions can reach 10p • 

All models of the formation of the quark-gluon plasma in nuclear collisions 
require the information about its rate of thermalization. Does the matter inside 
a fireball thermalize sufficiently fast, so that a thermodynamic description makes 
sense? For this, the thermalization time must be shorter than the fireball lifetime 
Tf ~ 10~ 22 sec . The thermalization time, or the time of local equilibration, Ti oc , 
can be estimated as follows [13-15]. The local equilibrium time writes as r ioc = A/c , 
in which A = (po") -1 is the mean free path of a particle in a medium; p = a -3 , 
the density of matter; a , average interparticle distance; a ~ b 2 , crossection; b , 
interaction radius. Accepting the values a ~ b ~ lfm typical of nuclear matter, 
we have A ~ lfm and T\ oc ~ lCT 23 sec . Because of the inequality T\ oc <C Tf , the 
thermalization inside a fireball is likely to be reached. 

The possibility of speaking about the thermodynamics of strong interactions, 
separately from electromagnetic and weak interactions, is based on the fact that it 
is just this type of interaction which in many cases plays the dominant role. Really, 
the dimensionless coupling constant of strong interactions a s ~ 1 is much larger 
than the coupling constants of electromagnetic interactions, a e as 1/137 ~ 10~ 2 , 
and of weak interactions, a w ~ 10~ 5 , to say nothing of the coupling constant of 
gravitational interactions, a g ~ 10~ 12 . The corresponding interaction times, at the 
energy 1 GeV characteristic of high-energy physics, are r s ~ 10~ 24 sec for strong 
interactions, r e ~ 10~ 21 sec for electromagnetic interactions, and r w ~ lCT 10 sec 
for weak interactions. Therefore, during the lifetime of a fireball Tf ~ lCT 22 sec 
electromagnetic and weak interactions do not play any role. In the interval of time 
T[ oc < t < Tf inside a fireball there may exist an equilibrium state of strongly 
interacting particles. 

The sole consistent way of calculating thermodynamic characteristics in the 
frame of quantum chromodynamics is perturbation theory, which is quite similar 
to that of quantum electrodynamics [16,17]. However, the effective coupling param- 
eter of strong interactions becomes small only at asymptotically high temperatures. 
In the most interesting region of temperatures around @ c ~ 160 MeV , where 
deconfinement occurs, the coupling parameter is large and perturbation theory does 
not work. For describing the whole range of thermodynamic variables several sta- 
tistical models have been suggested . 

2 Statistical Bootstrap Models 

The usefulness of applying statistical methods for considering the heated and com- 
pressed nuclear matter has been understood long time ago. Let us mention, e.g., Fermi 
[18]. 

The first statistical model of nuclear matter under extreme conditions, such as 
being realized inside fireballs, has been proposed by Hagedorn [19] (see also [20]) 
and called the statistical bootstrap model. In this approach it is assumed that, at 
zero baryon density ub = , various hadrons can be generated from vacuum with 
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the mass distribution 



p(m) = pdisim) + pconim); m G [0, oo), 

in which the first and the second terms correspond to discrete and to continuous 
mass spectra respectively, 

Pdisim) = J2dH m - m) [1 - 0(m - m )] , 



Pcon(m) = Q(m - m )-^ exp (J^j 



p = © J p(m) exp (^—f3Vk 2 + m 2 ) - d dm, 



where ©(•) is the unit-step function; Q , a degeneracy number for spin-isospin 
states; and the parameters are 

m = 1000 MeV, a = 6.5 x 10 3 MeV 3/2 , 6 = 160 MeV. 

The pressure for the ideal hadron gas is written in the classical Boltzmann approx- 
imation, 

d k 

where is temperature in energy units, and j3 is inverse temperature, /3Q = 1 . 

As can be easily checked, the pressure in this model diverges for all O > 0o . 
From here it was concluded [19] that Oo is the limiting temperature of the Universe. 
The concept of the existence of a maximal temperature of the Universe is, of course, 
quite artificial, therefore another interpretation of this divergency of pressure has 
been proposed [20] treating Oo as the deconfinement temperature. 

After Collins and Perry [21], the geometrical scenario of the deconfinement be- 
came popular [22] . According to this, the increase of temperature or baryon density 
leads to the rising number of hadrons. The latter are assumed to have finite volumes 
[23] . When the number of hadrons becomes so high that their close packing occurs, 
then they fuse into one gigantic hadron occupying the whole system. The Hage- 
dorn temperature Oo is interpreted as the fusion temperature, and the gigantic 
hadron cluster is identified with the system in the quark-gluon-plasma state. This 
geometric scenario reminds the percolation transition [24]. 

Following the geometric interpretation, the bootstrap model was modified [20] 
invoking the excluded-volume approximation. In this, one considers particles 
having the volumes v±, t>2, . . . , v n as moving in the free volume 

N 

v N = v-J2vj, 

where V is the total volume of the system. The pressure in the excluded-volume 
approximation becomes 



O 00 1 



V N J p{m) exp (—fiVk 2 + m 2 ) -^^—dm 



(27T)* 
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where again the Boltzmann approximation is also used. The mass distribution for 
discrete spectrum is taken in the same form as above, and for continuous spectrum 
it is slightly modified as 

Pcon(m) = 9(m - m )-^-exp ( ^ ) , 

with | < a < | . Now, the pressure is everywhere finite and positive becoming 
zero at the same temperature 9^ ~ 60 . The temperature 9^ , where p(Od) = 
is interpreted as the temperature of hadron fusion into a gigantic cluster. However, 
the thermodynamics of the system at 9 — > 00 has nothing to do with that of the 
ideal quark-gluon plasma. 

To overcome the latter deficiency of the model, it has been argued that taking 
into account hadron compression can save the situation. This can be done [25] by 
complicating the mass distribution writing its continuous part as 

p con (m,v) = 9(m - b v - m )Q(v - v )a (m - b v) a x 




which contains now seven fitting parameters: mo, do, a, 60, ^o, 7, °o • 

Now, in accordance with the geometrical scenario, the number of hadrons at low 
temperature 9 < 9<f is proportional to the free volume, and at high temperatures 
9 — > 00 this number tends to one symbolizing the formation of a gigantic cluster. 
However, at the 1-st order transition temperature 9^ the number of hadrons 
N(Qd, V) diverges for any finite volume V , which is unreasonable. 

The bootstrap models, in addition to the arbitrariness in postulating the mass 
distribution, contain internal deficiencies leading to the existence of instabilities 
contradicting the necessary stability conditions for statistical systems [26,27]. 

3 Pure Phase Models 

An evident idea would be to follow the standard Gibbs prescription for consider- 
ing phase transitions between two phases. Treating the deconfinement as such a 
phase transition, one assumes the existence of two types of pure phases. At low 
temperatures and baryon densities this is a pure hadron phase with features typical 
of the normal nuclear matter [28] , and at high temperature or baryon density this is 
the quark-gluon phase described by perturbative QCD [17]. Let Q± be the grand 
potential of the quark-gluon plasma, and Q2 , that of the hadron matter. Accord- 
ing to the Gibbs rule, a phase transition occurs when Qi(9,//b) = f2 2 (9,/i£) , 
where /i# is the baryon chemical potential. This equality gives a transition line 
9cz = 9rf(/iB) • Expressing here the baryon potential hb = /Ub(^_b) through the 
baryon density, we may write 9^ = ®d{ n B) ■ 

The possibility of using perturbation theory for the high-temperature quark- 
gluon plasma is based on the property of asymptotic freedom. According to this 
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property, the running coupling constant 



[11N C - 2N f ) \n(q/A) 



511nln(g/A) 
1211n(g/A) 



in which N c and Nf are the number of quark colours and flavours, respectively, 
A pa 200 MeV is a scale parameter, and q is momentum, tends to zero as q — > oo . 
In the integral over momenta, defining the grand potential, the main contribution, 
when — > oo , comes from q ph Q . Hence, it is possible to get an expansion in 
powers of a s (Q) = g 2 (0)/47r with the effective coupling parameter 



A®) 



24tt 2 



(llN c -2N f )\n(Q/A)' 



As a result of this expansion [29], neglecting quark masses, one has for the grand 
potential 

£ = + B, 

in which a nonperturbative term B is included and the notation 

A = A + A 2 g 2 + A 3 # 3 + A 4 # 4 In g 

is used, where 
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/i/ being the chemical potential of an / -flavour quark. Note that the term 0(g 6 ) 
cannot be calculated by perturbation theory because of unrenormalizable infrared 
divergences. 

Take into account that the number of quark colours is N c = 3 , and consider 
for simplicity the case of zero baryon density % = , so that /if — . Then the 
pressure of the quark-gluon plasma is 



Pi 



B, 



6 



with the expansion coefficients of A being 



45 V 32 V 2 6 V 12 f 1 ' 



9 / 1 \ 3 / 2 ^ 

A 3 = — 1 + -N,) , .4.= 



(i + ijv,). 



3tt V 6 V ' 2tt 2 
The low-temperature hadron phase is often modeled [4,7,30] by a gas of massless 
noninteracting pions, which for the grand potential yields 



The corresponding pressure is 



^2 = 

V ~ 30 ' 



^2 VT 2 . 

P2 = -y = 3O - 



d ( f2 \ dp 
s = ' 



Equating Qi with Q 2 , or pi with p 2 , one obtains the deconfinement 
temperature 

8, = ^; 7.(x-3" 1/4 . 

To check the phase transition order, let us find the latent heat at the transition 
temperature. 

Define the energy density 

e = s© - p + plbtib 

and the entropy density 

" <)(-) \V) d& 

The latent heat density is 

Ae d = £l -e 2 = Qd&s d (e = e d ), 

where As rf = Si — s 2 is the entropy density jump at = 0<f . 
For the considered case we have 

2tt 2 

s 1 = (4A + C)Q 3 , s 2 = —Q 3 , 

lo 

with the notation 

C- Q dA 33-2^ 
6= °90 - 24vr 2 5 

where the equation 



A 2 + ^A 3 g + ^A 4 (l + 4\ng)g 2 



dg 2 33 - 27Vj 4 
90 ~ 24vr 2 5 
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is taken into account. The energy densities are 
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ei = (3A + c)e 4 + B, ^ = ^e 4 



Thus, for the latent heat one gets 

Ae d = AB (l + ^ 7 4 

The nonperturbative term B is usually treated as the bag constant, with B 1 / 4 
ranging in the interval (150 300) MeV . For estimates, we may accept the value 
B l / A = 225 MeV from the middle of this interval. 

One often assumes that the quark-gluon plasma is an ideal gas of free quarks 
and gluons, that is g = . If this is so, then for Nf — 2 we have 7 = 0.72 . 
The deconfmement occurs at = 162 MeV being a first-order transition with 
the latent heat Ae d = AB 1 GeV / fm 3 . The found deconfmement temperature 
almost coincides with the characteristic energy E = 157 MeV discussed in 
Introduction. This picture would seem quite reasonable if it would not be absolutely 
wrong. Really, the effective coupling <?(©) — > 00 for temperatures close to the 
scale parameter A m 200 MeV . Therefore, the assumption that, in the vicinity 
of the deconfmement point, the quark-gluon plasma is an ideal gas of quarks and 
gluons is senseless. This would happen only at temperatures at which g 2 (0) <C 1 , 
that is when 

e>>Aexp l^V/J- 

The latter inequality becomes valid only at very high temperatures 6 3> 10 6 MeV . 

In this way, nonperturbative effects around the deconfmement transition are very 
strong, and it is not correct to try taking account of them by the simple addition to 
the grand potential of a term B . 

Some nonperturbative effects can be included into consideration by resorting to 
the effective-spectrum approximation [31,32], when one postulates for the spectra 

of quarks and gluons the form Si{k) = \Jk 2 + mf + Ui , in which k is the modulus 
of momentum; , a mass; C/j , an effective mean field; and % — q,g enumerates 
quarks and gluons. This approximation yields the results similar to the ideal gas 
picture: The deconfmement is a first-order transition occurring at 6 ~ 160 MeV . 

The reason why the effective-spectrum and ideal-gas approximations are close 
to each other can be understood as follows. The effective-spectrum approximation 
may be interpreted as a result of a renormalization of perturbative series. As an 
example, we may use the self-similar renormalization [33-36] differing from other 
resummation techniques by the possibility of checking its range of applicability at 
each step. Consider the coefficient A in the grand potential Qi of the quark-gluon 
plasma as an effective limit of the sequence {_/&(<?)} with the initial approximation 
fo(g) = A + A 2 g 2 , the first approximation fi(g) = fo(g) + A 3 g 3 , and so on. The 
simplest variant of the self-similar renormalization [33-36] gives the renormalized 
coefficient 

,*_ a 4A l9 2 

° + (A 3 g-2A 2 y 
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This quantity, as A 2 < , is finite for all g including g — > oo . In the whole 
diapason of g G [0, oo) the deconfinement temperature does not change much: 
for g — , with Nf — 2 , it is ©d = 162 MeV , while for t/ — > oo it is 
0^ = 176 MeV . However, the renormalized value A* is truthful only when the 
self-similar renormalization is stable [37-40]. For this we need that the corresponding 
mapping multiplier Mi(g) and the Lyapunov exponent Ai(g) would satisfy the 
stability conditions: \Mi(g)\ < 1 and Ai(g) < . In the considered case 

OA 

M 1 (g) = l + A 1 (g), A^g) = ^g. 

The Lyapunov exponent, since A 2 < and A 3 > , is always negative. The 
condition M 1 (g) holds only for g < 4\A 2 \/3A 3 , that is for g < 1 . 

Thus, the renormalization of the grand potential, starting from the perturbative 
expression, does not essentially change the results. And it is clear why: Really, the 
ideal-gas picture, which makes the basis of the perturbative expansion, contains 
no information on bound states that should appear as the coupling constant g 
increases. The situation with the quark-gluon plasma is somewhat similar to that 
of the electron-ion plasma in which there can exist bound as well as free electron 
states [41]. 

The value of the deconfinement temperature obtained in pure-phase models can 
be quite reasonable, in the same way as the simple estimate of Introduction is 
such. However, it would be hard to believe that these models can correctly describe 
the character of the deconfinement transition and the behaviour of thermodynamic 
functions. 



4 Lattice Numerical Simulation 

The idea of using a discrete space-time lattice to regularize quantum field theories 
opened the entire repertory of statistical physics for the analysis of nonperturbative 
properties of these theories. The application of Monte Carlo simulation techniques 
turned out to be a powerful approach allowing to perform a quantitative study of 
nonperturbative aspects of quantum chromodynamics. The lattice reformulation of 
QCD has been described in several surveys (e.g. [4,29,42]), therefore below we only 
slightly touch the principal points of this approach. We shall mainly discuss the 
predictions of the lattice QCD and its simplified variants for the deconfinement 
transition. 

The first step towards the finite temperature study of QCD consists of intro- 
ducing the imaginary time t = —ir and of rewriting the partition function as a 
path integral of the exponential of the Euclidean Lagrangian density over all fields 
in the problem. The second step defines the cubic four-dimensional lattice with the 
sites x = {x,t} where x£ Z 3 is a real-space lattice vector and r e Zi is 
conventionally called the temporal variable. The lattice spacings in the spatial and 
temporal directions are denoted by a a and a T , respectively. If N a and A^ T are 
the number of sites in the corresponding directions then the volume V and the 
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temperature 9 of the system are given by V = (N a a a ) 3 and (5 = N T a T , where 

p = e- 1 . 

A gauge-invariant theory on the lattice is usually formulated in terms of link 
variables U£ and site variables ip(x) and ip(x) . The link variable U£ is 
associated with a link leaving site a; in a direction fi = 1, 2, 3, 4 and it is a matrix 
U£ G SU (JV C ) in the space of colour indices, N c being the number of colours. The 
site variables ip(x) and ip(x) are associated with each site x of the lattice and 
they carry colour, flavour and spin indices. Also, ip(x) and ip(x) , representing 
fermion fields, are treated as Grassman variables. The link and the site variables 
satisfy the periodicity conditions 

Ut fl = U^, i/;(x,0) = -rl>(z,j3). 

Then, in terms of these variables, the partition function is written in the form 
of a path integral 

z = /nnn e ^w(^)#(^), 

J X fl x > 

in which the action 

S = S(U^^(x')^(x')) = S G + S F 

consists of a gauge action Sq and a fermionic action Sp • The partition function 
and other thermodynamic functions are calculated by using the Monte Carlo pro- 
cedure on a finite lattice of N% x iV r sites, the maximal number of sites in each 
direction being around N a = 24 and iV r = 24 . 

The deconfinement transition has been studied, in the frame of the lattice QCD 
or its pure gauge variants, by many authors. Here we cite only some review-type pa- 
pers [43-52]. Note, that, in addition to the confinement transition, another transition 
related to the spontaneous breaking of chiral symmetry can occur. The mechanisms 
leading to these transitions are seemingly unrelated and it thus has been speculated 
that QCD may undergo two separate phase transitions. However, the Monte 
Carlo data for pure gauge models and for the variants with massless quarks suggest 
that these two transitions coincide. And if the finite masses of quarks are taken 
into account then the chiral symmetry, strictly speaking, happens only in the limit 
— > oo . In what follows we shall speak solely about the deconfinement transition. 
One of the main discoveries of the lattice simulation has been the fact that the 
deconfinement can be of quite different character for different systems. 

SU(2) Pure Gauge Model 

The common convention is that the deconfinement for such quarkless models is a 
second-order transition occurring at ~ 210 MeV . The errors in calculating the 
pressure and energy density, related to the finiteness of the used lattice, are around 
5% above and 30% below Q d . 

SU(3) Pure Gauge Model 
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The deconfinement has been found to be a first-order transition at 0^ pa 
225 MeV with the latent heat Ae d ~ \£i(Qd) • The errors of calculating the 
pressure, energy density and entropy are about 10% above and 30% below Q d . 

SU(3) Model with Quarks 

In the presence of dynamical fermions the situation is by far more complicated 
as it is difficult to perform high statistics analysis of full QCD on large lattices 
that would allow a detailed finite size study as it has now been done in the pure 
gauge sector. Monte Carlo simulations for full QCD with dynamical quarks have 
by now been performed only for the case of zero baryon density. The transition has 
been found to be rather sensitive to the choice of quark masses and the number of 
flavours. It seems that for Nf > 4 the transition is first order for all quark masses. 
For Nf < 4 the situation is still to some extent uncertain. Nonetheless, there are 
strong indications that for physical quark masses and Nf = 3 the transition is 
likely to be a continuous crossover, occurring at Q d m 150 MeV . However, one 
should once again stress the uncertainties in the determination of the order of the 
transition in Monte Carlo simulations, with quarks, even at zero baryon density. 

What has been found in the lattice simulations with certainty [52,53] is that non- 
perturbative effects are quite strong around the deconfinement transition persisting 
till about 2<d d . Also, as is seen, no statistical bootstrap or pure phase model is 
able to describe the variety of different transition orders discovered in the lattice 
simulations, since these models always predict a sharp first-order transition. 

5 Dynamical Confinement Scenario 

Correlation functions and related susceptibilities are one of the key tools used for in- 
vestigating phase transitions [54]. This concerns as well QCD correlation functions 
[55]. The latter have been intensively studied in lattice numerical simulations [56- 
60]. The results of these studies revealed nontrivial effects in the high temperature 
phase connected with strong quark-antiquark correlations, which was interpreted 
[56-60] as the existence of hadronic modes, even at high > Q d ■ De Tar [56] 
proposed that the high temperature phase might be dynamically confined, in the 
sense that the long range fluctuations are colour singlet modes, and that the poles 
and cuts in the linear response functions of the hadronic phase go over smoothly 
into those of the high temperature one. This scenario of dynamical confinement 
presupposes that, actually, there is no transition at all, and that there is only a 
smooth crossover between hadronic matter and the quark-gluon plasma. 

It is conjectured [56] that the characterization of the plasma as a weakly in- 
teracting gas of quarks and gluons is valid only for short distances and short time 
scales of the order 1/0 , but that at the scales larger than l/g 2 Q , where g 2 
is the running QCD coupling, the plasma exhibits confining features similar to 
that of the low-temperature hadronic phase. The confinement scale, for instance, 
at ~ 200 MeV , is expected to be of order 1 fm and 10 _24 sec . 

Hatsuda and Kunihiro [61-64] considered the density-density correlation func- 
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tions of quarks for the Nambu-Jona-Lasinio model with a QCD -motivated ef- 
fective Lagrangian. They found precursory collective excitations existing in the 
high-temperature phase and corresponding to correlated qq pairs. This means 
that meson modes do exist above as well as below transition temperature. Really, 
all poles and cuts of the correlation functions in the momentum-energy representa- 
tion are the same at all temperatures. Thus, quarks, antiquarks, and gluons should 
also exist in the low-temperature as well as in high-temperature phase. Such pre- 
cursor effects are analogous to pretransitional fluctuations in superconductors [65] 
or superfluid He 3 [66]. 

The gradual change of the excitation spectrum from hadronic states to quarks 
and gluons, and the survival of hadronic modes above the transition temperature, 
have been confirmed for QCD in the instanton-liquid approach [67] and in the 
magnetic-current approximation [68]. In the latter case, quarks are assumed to 
interact at high temperature solely through magnetic current-current interactions, 
the electric ones being screened. 

The following picture [68] may serve as an intuitive illustration of the dynamical 
confinement [56-60]: The current-current interactions persist above 0^ and force 
any quark-antiquark pair (and, may be, every three-quark state) to correlate into 
colour singlets. As the quarks are moving in the heat bath, the string connecting 
them for colour neutrality are constantly breaking and reforming, which can be 
interpreted as hadrons going in and out of the heat bath. 

If the transition between the hadron matter and the quark-gluon plasma is a 
gradual crossover, as follows from the dynamical confinement scenario, then why in 
many lattice simulations this transition is found to be of first order? The answer 
to this question was suggested by Kogut et al. [69] explaining the first order of this 
transition merely as due to the finite size of the lattice. Calculating thermodynamic 
characteristics, such as energy densities, and correlation functions for several lattices 
with the sizes 8 3 x 4, 12 3 x 4, 12 3 x 6, and 16 3 x 4 , it has been found [69] that the 
quark-gluon plasma transition becomes less abrupt as the lattice size increases and 
the evidence for a first order transition becomes weaker. Lengthy runs for N a = 16 
showed no evidence for metastability, so that Kogut et al. [69] concluded that their 
results suggest that there is no abrupt transition at all, but only a smooth crossover 
phenomenon. 

To describe the coexistence of gluons and glueballs above for a quarkless 
SU(2) system, De Grand and De Tar [57,58] proposed to write the grand potential 
as the sum 



o 

where Q\ corresponds to gluons and Q2 , to glueballs; Q is a degeneracy 
factor; the Bose distribution rii(k) = [exp{/?£j(A;)} — 1] _1 contains the spectrum 
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£ i(k)= Jk 2 + mf with the mass m,\ = for gluons and the mass m 2 1000 MeV 



for glueballs. The cut-off momentum k is postulated, to fit lattice data, to be 
k -> oo below d and A; = 2.869^/(0 - 9 d )] a3 above B d . This 
phenomenological model describes only the second-order transition in pure 577(2) 
systems. The temperature dependence of the cut-off momentum A; is too arbitrary 
to permit a straightforward generalization to systems containing quarks as well as 
various hadrons. 

The main value of the dynamical confinement scenario is the clear understanding 
that hadron and quark-gluon degrees of freedom, generally, coexist. 



A model in which quarks could coexist with three-quark clusters, that is, nucleons 
has been suggested by Clark et al. [70] for zero temperature and considered by Bi and 
Shi [71] at finite temperatures. Nucleons and quarks in this model are intermixed 
inside the same system. This kind of mixture should be distinguished from the Gibbs 
mixture, in which different phases are separated in space having only a common 
interphase boundary (see [15]). The quark-nucleon mixture is rather similar to a 
binary liquid mixture [72]. A two-component system can, generally, stratify in space 
if there are no chemical reactions between components. The quark-nucleon mixture 
does not stratify because of the possibility of formation and decay of nucleon clusters, 
which is taken into account by the relation pn — 3/i ? between the nucleon, p^ , 
and quark, p q , chemical potentials. Since it is assumed that inside a nucleon bag 
and outside it the QCD vacuum is different, the nonstratified mixture corresponds 
to nonuniform, or unhomogeneous, vacuum. In some sense it reminds the twinking 
vacuum discussed with respect to the instanton-liquid picture of the QCD vacuum 
[29]. A nucleon can be interpreted as a droplet of a denser phase inside the rarified 
plasma phase. Nucleons among quarks are alike fog drops in air. Another analogy 
is magnetic bubbles in magnets [73]. 

Clark et al. [70] used the excluded-volume approximation when particles move 
not in the whole volume V but only in the free volume V — V — N n Vn , where 
Nn is the number of nucleons of the volume vn each, defined by the bag model. 
This renormalizes the volume V — > t^V by the factor £ = Vq/V = 1 — pnVn , 
where pn = N^/V is the nucleon density. Then, e.g., the baryon density takes the 
form 



in which, in addition to the geometrical interaction taken into consideration by 
the factor £ , the perturbative quark-quark interaction is simulated by the factor 
1 — 2a s /ii . The Fermi distributions rii(k) contain the free spectra u)i(k) = 



6 Clustered Quark Model 





+ mf — pi . The energy density reads 



s = e q + En + B, 
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e q = eC 9 (l - / e q (k)n g {k)j^- 3 , 

6 N = ^Cn J £N(k)n N (k)-^J^, 

where the nonperturbative energy is given by a bag constant B . 

Comparing the grand potentials for this model with B l / A = 171 MeV and those 
of pure quark and of pure nucleon phases, it was found that the mixed clustered 
matter is more profitable at n B > 8n 0B and < 50 MeV . Of course, one 
should not take too seriously the numerical predictions of this model which is too 
oversimplified to be realistic. For instance, pions and gluons are not considered here 
although they should play an important role in the deconfinement transition. Thus, 
pion fields make a nucleon bag unstable at quite moderate density and temperature 
[74,75] close to those discussed in Introduction. 

The most valuable in the clustered quark model is the idea of nonuniform vac- 
uum containing such local fluctuations that nucleons can be formed floating in the 
surrounding quark matter. 



7 String-Flip Potential Model 

An extreme realization of the dynamical confinement scenario is given by the string- 
flip potential model [76-81] based on the quantum-mechanical Hamiltonian 

H(r 1 ,...,r N )=j:i [ -^-j + V(r 1 ,...,r N ) 

of nonrelativistic quarks. In the interaction term 

V(r 1 , . . . , r N ) = min ^ $(r^) (r^ = | r t - 7j |) 

<ij> 

the minimum is taken over all ways to group iV quarks into hadrons, so that 
the summation actually goes only over < i,j > pertaining to the same hadron 
consisting of 2, 3, 6, 9, 12 or 15 quarks. Using the variational wave function 

*(r i, . . . , r N ) = exp{-/3V(r i, . . . , r N )}D(r l , . . . , r N ), 

in which (3 is a variational parameter and D(. . .) is a free Fermi gas Slater 
determinant, one minimizes the energy E((3) = H^>)/(^>, ^) with respect to 
(5 . This gives (5 — (3{p) as a function of density. 

Variational Monte Carlo calculations for this model have been accomplished with 
the harmonic, $(r) = |mu; 2 r 2 , and linear, $(r) = or , confining potentials. The 
wave function has been symmetrized with respect to all of the quark coordinates, 
including these pertaining to different hadrons. This takes into consideration the 
exchange quark energy which is found [82] to be about 100 MeV . 
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The results show that the quarks always coalesce into the lowest energy set of flux 
tubes, which is characteristic of an adiabatic approximation to the strong coupling 
limit of QCD . At low densities p < p quarks cluster into isolated hadrons. As 
the density increases, the value of (3(p) decreases, first slowly, but at p 1.5po 
exponentially. For p > 2p Q the wave functions of separate hadrons strongly overlap, 
which may be interpreted as a transition to the quark matter, although formally all 
quarks are yet confined. The complete deconfinement occurs at asymptotically high 
densities, when (3(p) — > , and the wave function becomes that of a free Fermi gas 
of quarks. 

The quark-quark pair correlation function was calculated [76,79] using Monte 
Carlo techniques and the longitudinal response function, using molecular dynamics 
simulations [83]. At low densities these functions have the properties typical of a 
set of isolated hadrons. For instance, the pair correlation function has a sharp peak 
displaying strong correlations between quarks. The form of this function is similar 
to that of liquids [84,85]. As the density increases, the peak in the pair correlation 
function relaxes smoothly to that of a Fermi gas of quarks, disappearing completely 
as p — > oo . 

Thus, the string-flip model demonstrates a smooth transition from hadron mat- 
ter at low densities to a free Fermi gas of quarks at high densities. Quarks at all 
densities are to some extent confined becoming an absolutely free gas only in the 
limit p — > oo . 

The main difficulty in dealing with this model is the necessity of the minimiza- 
tion of the Hamiltonian interaction term over all strings connecting quarks, which 
involves complicated Monte Carlo calculations. This procedure is necessary for elim- 
inating colour van der Waals forces between colour singlet hadrons. These forces 
are present in confining two-body potentials and, if not eliminated, produce large 
spurious energies in nuclear matter. 

However, in a consistent statical approach this kind of divergence can be auto- 
matically cancelled if one takes into account the corresponding correlation functions, 
called smoothing or screening functions, which smooth the interaction potentials in 
the regions of their divergence [86-88] . It is possible, starting with divergent inter- 
action potentials, to construct regular iterative theory [35,89] for Green functions, 
whose each step involves only smoothed potentials containing no divergencies. In 
the lower orders of this correlated interaction theory [35] the situation is such that 
could be obtained by replacing from the beginning the bare interaction potential 
by a smoothed one, that is by defining an effective Hamiltonian with a pseudopo- 
tential instead of the initial potential [86]. Such a replacement is an approximation 
neglecting some double and all triple correlations [35]. 

The Hamiltonian of the string-flip model including the minimization over strings 
is an effective Hamiltonian, in which this minimization plays the role of smoothing. 
In lieu of this complicated minimization requiring Monte Carlo techniques, it is 
equival ent to replace the confining potential by a smoothed one. 

Ropke at al. [90-93] have developed a many-body approach to quark-nuclear 
matter generalizing and simplifying the string-flip model. This approach starts with 
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an effective Hamiltonian 

N ( V 2 \ 1 N 
H eff (r r N ) = E k" + 3 E ^//(^ij), 

in which stands for the quark masses, and the effective interaction V e ff(r) = 
s(r)$(r) is screened by a function s(r) defining the probability that two quarks 
at distance r are next neighbours. The screening function is to be found from 
some additional equations. In the case of independent particles with density p , 
one has s(r) = exp(— 4^pr 3 ) . Inside a hadron, when r <^ a , where a = p -1 / 3 is 
an average interquark distance in the system, we have s(r) ~ 1 , so quarks interact 
through a bare confining potential. For the quarks pertaining to different hadrons, 
i.e. when r > a , we get s(r) ~ , which is called the saturation property [91]. 

The many-body approach to clustering quark matter [90-93] makes it possible 
to give a unique description of both the clustered hadrons as well as the free-quark 
phase in the same way as it has been done for the nuclear matter whose nucleons 
could form the deuteron, triton , and a -particle bound states [94,95]. The total 
density of clustering matter is a sum p = J2i z iPi , i n which Zi is a compositeness 
factor and p is the density of an i -component of clusters composed of Zj particles 
each. 

In the many-body approach the quark system may contain both quasifree quarks 
(scattering states) and various clusters (bound states). At zero temperature, nucle- 
ons start dissolving at 113 ~ 4.3nos • At zero baryon density the transition of 
hadron matter to free-quark phase starts at O^ ~ 200 MeV where some of meson 
states are already dissolved. But other bound quark states do exist yet at this tem- 
perature. The dissolution of the bound states occurs gradually with temperature. 
Some mesons survive up to ps 20 d . The gradual dissociation of bound states 
above 0^ is in agreement with the dynamical deconfinement scenario. 

Although the many-body approach provides the principal possibility of treating 
bound states together with quasifree ones, it has the following technical obstacles. 
Each bound state is described by a Bethe-Goldstone-type equation for a two-or 
three-particle Green function, depending on the compositeness of this bound state. 
This type of equations, as is known, is very difficult to deal with. When there are 
many bound states interacting with each other as well as with unbound particles, 
then one has to deal with a system of many interrelated Bethe-Goldstone equations, 
in addition to several equations for one-particle Green functions. In such a case the 
problem becomes practically as difficult as chromodynamics itself. 

8 Concept of Cluster Coexistence 

As follows from the previous Sections, there is need of such a statistical model that 
could provide a realistic description of clustering matter being at the same time 
treatable. One almost evident simplification would be, instead of suffering with a 
system of Bethe-Goldstone equations, to consider bound clusters as separate objects. 
This can be understood as a kind of renormalization by integrating out the internal 
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degrees of freedom corresponding to the motion of particles inside each of clusters, 
so that only the center-of-mass degrees of freedom are left. 

Enumerate all possible cluster states, including those of unbound, or free, parti- 
cles, by the index % . Let ipi( r ) = [ l l J ?{ r )} be the field operator of an % -cluster; 
this operator being a column in the space of quantum degrees of freedom indexed 
by a . The quantum degrees of freedom include spin, isospin, colour and flavour 
indices. The field operators satisfy the commutation relations. 

mr),4>J(r)] T = 0, [1>?(r)J] (7% = 5^5(7 - 7), 

where the upper sign stands for Bose and the lower, for Fermi statistics. 
Take an effective Hamiltonian in the form 

H = E-J2^N l , (1) 

i 

in which the energy operator 

£ = E^ + E4- (2) 

i i,j 

consists of the single-cluster energies 

E t = f4{7)K i {V)A{r)d7, (3) 

where .fQ(V) is a kinetic-energy operator, and of the interaction energies 

= ^J ^1(7)^(7)^(7 - 7')^ j (7')i) i (7)d 7 d 7', (4) 

with the interaction potentials $ij(r) having the symmetry property 

Sy(f) = <M- 7) = (5) 

In the second term of (1), /i, is the chemical potential of % -clusters and 

Ni = J 4{r)i>i(7)d7 (6) 

is a number-of-cluster operator. 
The grand potential 

n = -einTr e- pH (7) 

defines the pressure 

P=~ = finite-^. (8) 
The energy density is given by 



e = ^{$), (9) 
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where (...) denoted the statistical averaging. 

For each kind of clusters we may define their density 



the number of clusters 



and the average density 



Pi(?) = (4(rM(r)>, (10) 
N i = (N i ) = j Pi (r)d7, (11) 

Pi = ^ = - j Pi {r)dr . (12) 

The number of constituents of an % -cluster is called the compositeness factor 
Zi . The total number of elementary particles is 

N = 5>JVi, (13) 

i 

so that the average density of the system can be written as 

TV 

P= y = Y, Z iPi- ( 14 ) 

i 

The chemical potentials Pi can be defined from the conservation laws accepted 
for the given system. For example, if the numbers iVj for each sort of clusters are 
fixed, then (12) gives P i(Q, P i) . If the total number of constituents (13) is fixed, 
then from the equilibrium condition SQ = we have 

to = Mi > (15) 

Z\ Zj 

which, together with (14), defines Pi (Q )P ) . When neither N { nor N are fixed, 
then ^ = . 

An important quantity is the cluster probability 

Wi = z i— = Zi— (16) 

enjoying the conditions 

o<Wi<i, E^ = 1 ? ( 17 ) 

i 

following from (14). With (12), we may rewrite (16) as 

Zl J Pl (7)d7. (18) 



Wi - 

N 



If the system of particles can form a number of bound states, then how could 
we define numerous interactions between clusters? The number of such interactions 
can be drastically reduced for those clusters whose interaction potentials $ij(r) 
have the same type of behaviour at large distance, for example, diminishing as r 
increases. Consider such clusters with a similar behaviour of ®ij(r) . Let two 
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clusters, m and n , coalesce into one, % , so that the compositeness numbers 
involved in the reaction m + n — > i satisfy the obvious relation z m + z n = Zi . 
And let aside be another cluster-spectator, j , as is shown in Fig.l. Assume that 
the coalescence of the cluster-actors does not influence the cluster-spectator, in the 
sense that its interaction with the initial two clusters is the same as with that one 
formed after the coalescence: 

* mj (r)+* nj (r) = *i3(r). (19) 
From this assumption, using the conservation law z m + z n = Zi , we find the relation 

®ij ( r ) _ ®mn( r ) 

ZiZj z rn z n 

permitting to express the interaction potentials of different clusters through one 
calibration potential. 



9 Screening of Interaction Potentials 

The interaction potential <&jj(r) may be divergent for some regions of r . However, 
it would not be correct to replace this potential by a smoothed, or screened, one 
just in the Hamiltonian (1). Such a screening appears in the process of decoupling 
of binary propagators with taking account of interparticle correlations, as is done in 
the correlated iteration theory [35]. If one wishes to write an effective Hamiltonian 
corresponding to a given correlated approximation, it is not sufficient to merely 
rearrange the operator terms of the Hamiltonian, but a nonoperator correcting term 
must be added. 

Let us illustrate this for the correlated mean-field approximation [35] when the 
Hamiltonian (1) may be presented in the form 



H = J2H l + CV, 



H i = J $(r 
in which the correlated mean field 



i>i(r)d?, (21) 



Ui(r) =£ / iy (r - r) Pj (r)d 7' (22) 

3 J 

contains the screened potential 

*H (?) = Sijir^ai?) (23) 
smoothed by the smoothing function having the symmetry property 

s i j(r) = s ji (r)=s ij (-r). (24) 
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The nonoperator correcting term in (21) is CV , and this cannot be put zero, as 
will be shown below. 

Since the exact Hamiltonian (1) does not depend on cluster densities, varying 
the grand potential (7) with respect to Pi(r) we have 

I 5H \ 

= 0. (25) 



5pi(r) \5pi(r) i 

In order that the exact Hamiltonian (1) would be correctly represented by the ap- 
proximate Hamiltonian (21) requires that the latter must satisfy (25), which yields 

SC + hx(^k) = 0. (26) 



$Pi( r ) v j \$Pi( r ), 



Substituting here 



transforms (26) into 



5C 1 _ fSUj(r) ,^,-/ n 

2 / -^Pj(r )dr = 0. (27) 



Spi(r) V ■ J Spi(r) 
If the smoothing function Sy(r) does not depend on Pi(r) , then (22) gives 

5UA7 



8 Pj {r ) 

And Eq.(27) becomes 

5C 



$ij(r - r ). 



^ ~ r')Pj(r')d r'= 0. 



Spi(r) V J 

The solution of the latter variational equation, up to a constant that can be omitted, 
is 

C = -^77 £ / ~ r')Pi{r)pj{r')d 7 d 7' . 

Note that the Hamiltonian (21) with the obtained correcting term could be derived 
from (1) with the substitution 

^j(?)^j(?')^(^Vi(7) -> 

-> Sij (r - 7) {$(r)ip i (r)p j (r) + pi(7)ip} (7)^(7) - Pi(r) Pj (r)} , 

corresponding to the correlated Hartree approximation [35]. 

In general, the smoothing function Sjj(r) is dependent on Pi(r) . There- 
fore, the correcting term is to be found from Eq.(26) or (27). In any case, the 
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correcting term depends on thermodynamic variables through the densities Pi(r) . 
Thus, neglecting this term would disfigure the correct statistical description, and 
the behaviour of thermodynamic functions could be completely spoiled. 

Let us pass to a uniform system when Pi(r) — p^ . In principle, a multicompo- 
nent system can display a variety of nonuniform states related to the solidification 
of one or several components. For example, an ensemble of fully ionized nuclei can 
form a crystalline lattice immersed in a uniform electron background, which models 
the high-density matter of white dwarfs [96]. It may be that some heavy cluster 
components crystallize while others are liquid-like as it happens in superionic con- 
ductors [97]. It also may be that heavy clusters form an amorphous solid while light 
ones move in a conduction band as in glassy metals [98]. In the cores of neutron 
stars a Gibbs mixture can exist when in some volumes of space a lattice structure 
appears while others are filled by a liquid-like nuclear matter [99]. We leave aside 
all these possibilities considering in what follows only uniform systems. 

In the uniform case, the mean field (22) becomes 

E/i(r) = = E/i (28) 

3 

with the interaction integral 

(29) 




Instead of (25), we have 



The correcting equation (27), defining the correcting term, changes to 



(30) 



5C ^-^ SU-j 

^ + ?^ ft = ' (31) 

The field operators, for a uniform system, can be expanded in plane waves, 

v V k L J 

Then an % -cluster Hamiltonian in (21) transforms to 

H i = ^u i (k)at(k)a i (k), (32) 
k 

with the effective spectrum 

Ui(k) = €$) - p h = Ki(k) + U t . (33) 

The number of clusters (11) reads 

Ni = £< a lUK(£0) = Ci £"*£)> (34) 
k k 
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where d is a degeneracy factor, i.e. the number of quantum states, and 



rii(k) = \ exp 



/3uji(k] 



(35) 



is the Bose (upper sign) or Fermi (lower sign) momentum distribution. The cluster 
density (12) is 

d k 



Pi = G J ni(k) 
For the grand potential (7) one gets 



(2tt) 3 ' 



(36) 



fti = TQVQ / In 



The pressure (8) writes 



l±ni{k) 



P = J2Pi ~ C i 



d k 



(37) 



Vl = ±Q( t / In 



l±rii(k) 



d k 
(2^)3 " 



(38) 



The energy density (9) is 



£i = C J £i(k)rii(k] 
The cluster probability (18) becomes 



d k 
(27p' 



(39) 



(40) 



An additional simplification comes for isotropic systems, which is usually as- 
sumed, when Si{k) = Siik) , where k = | k \ ■ Then u>i(k) = ojiik) and 
ni{k) = Ui{k) . If the spectrum u>i(k) is such that the asymptotic properties 

A; 3 ln^(A;)^0 (k — > 0), 

u>i(k) — > oo (A; — > oo) 
hold true, then the cluster pressure and the cluster energy density reduce to 



Pi 



6tt 2 



k%(k)rii(k)dk, 



k 2 ei(k)rii(k)dk, 



1 2tt 2 Jo 
where e'^k) = dEi(k)/dk . 

These are the basic formulas which we shall use in what follows. 
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10 Coexisting Multiquark Clusters 



The concept of cluster coexistence has been applied first to nuclear matter consisting 
of different multiquark clusters. The interest to this problem was motivated by 
the A.M.Baldin commulative effect [1] and the related discussion of the possible 
existence of multiquark clusters in nuclei [100,101]. 

Since at high density or temperature relativistic effects play an important role, 
the kinetic term of the cluster spectra is taken in the relativistic form Ki(k) = 
\Jk 2 + mf , where m 8 is the cluster mass. A system of 3 -, 6 -, 9 -, and 
12 -quark clusters has been considered in the excluded-volume approximation [102- 
106]. The 3 - and 9 -quarks are Fermions, and the 6 - and 12 -quarks are 
Bosons. The Bosons with the lowest mass, that is the 6 -quarks, can drop down 
in the Bose-Einstein condensate, when u; 6 (0) = , which fixes /i 6 . The chemical 
potentials always satisfy (15). 

In the excluded-volume approximation the interaction between clusters is con- 
sidered geometrically by putting $jj zero, but replacing the total volume V by 
the free volume V , 

V^V = V-'£N i v i , (41) 

i 

where V{ are cluster volumes. This reduces the volume 

V^iV; £=^ = 1-5^, (42) 

i 

by a factor £ e [0, 1] . Equivalently, this can be interpreted as a reduction of the 
degeneracy factor 

Ci^Ci = ZSi- (43) 
The reduction factor in (42) and (43) can also be written as 



d k 



(2vr) 3 ' 

Thus, for the density of clusters, pressure, and energy density one has 

P = ±e£Ci/Mi±*(*)]^p 

e = ECi/ei(*W*)(^s- (44) 
The volumes of clusters are supposed to be related by the equation 

(45) 



vrii rrij 
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This permits to express all cluster volumes Vi = mjf 3 /m 3 through the 3 -quark 
volume t> 3 = 47rr|/3 with the radius r 3 = 0.4 fm of the nucleon core. 

The used multiquark parameters are given in Table 1. The 3 -quark is a 
nucleon. The 6 -quark mass m 6 = 1944 MeV corresponds to an average value 
over the masses of several light narrow dibaryons that are claimed to be observed 
in experiments [107]. The mass m 6 = 2163 MeV is taken from the bag-model 
calculation of Jaffe [108] and the 9 - and 12 -quark parameters are elicited from 
the bag model of Matveev and Sorba [109] . 

The 6 -quark probability depends on the value of the mass m§ as is illustrated 
in Fig. 2. The results for other cluster probabilities are displayed in Figs. 3 and 4, 
where m 6 = 2163 MeV . The probabilities of heavy clusters are always very small: 
w 9 < 0.1 , w 12 < 0.01 . 

As far as we limited here by the length of this review, we shall not discuss in 
detail the results of our calculations, which can be found in the cited papers. We 
think that the presented figures speak for themselves: it is better to see once than 
to listen hundred times. 

The coexistence of nucleons with 6 -quark clusters has been considered as well 
in the mean-field approximation [110] with the effective interaction potential 

= 2n^i-5(j) — — exp (—m n r) . 

The first term here is the Fermi pseudopotential for the core interaction with the 
scattering length = |(aj + %) , where = an , and the reduced mass 
rriij = miirij/ {m,i + rrij) . The second term is caused by the one-pion exchange; 
a n = 0.08 being the pion coupling parameter and m w = 140 MeV , the pion mass. 
Again, to reduce the number of model parameters, the relation af/m, = a^/rrij , 
similar to (45), is accepted. Then, all scattering lengths = a 3 (mj/m 3 ) 1 / 3 are 
expressed through the nucleon scattering length a 3 = 1.6 fm . The 6q -probability 
is shown in Fig. 5. 

Note that the mean-field approximation is valid if |C/j| <C rrii , which is true for 
the densities up to about 10p • 

The models of this Section serve rather as a qualitative illustration of cluster 
coexistence. They can have sense only at temperatures and densities much lower 
than those characteristic of deconfinement, as unbound quarks are not included here. 

11 Baryon Rich Matter 

Include into consideration unbound (quasifree) quarks that can, in principle, coexist 
with multiquark clusters. Whether and when quasifree quarks really coexist should 
be determined in a self-consistent way from the conditions of thermodynamic ad- 
vantageousness and stability. To compare the results with those of the previous 
section, consider again the case of baryon rich matter, when p = 3hb , that is when 
the generation of particles from vacuum can be neglected. 

Denote the chemical potential of a quark by /i = fi q . Then relation (15) yields 
Hi = Zifx . 



24 



The strengths of characteristic interactions between baryons and between quarks 
are of the order of or higher than the expected deconfinement temperature, so these 
interactions must be taken into account. The mean field acting on quarks may be 
written in a bag-model-motivated form [32] as U q = BV/N = B/p . This gives for 
quarks the spectrum 

£q (k) = yjk 2 + m 2 q + y (46) 

The mean-field term in (46) contains the total quark density p , which means that 
each quasifree quark interacts in the same way with other unbound quarks as well 
as with quarks entering into bound clusters. The interaction potentials between 
different baryons can be expressed, basing on relation (20), through the nucleon- 
nucleon interaction potential $33(7") , 

*v(r) = ^<I>„, (47) 

There are several such effective potentials obtained from nucleon-nucleon scattering 
experiments [111] or from analysing the deuteron properties [112]. We opt for the 
Bonn potential [113]. The common consensus is that thermodynamics of nuclear 
matter does not depend on the mutual orientation of spins of interacting nucleons. 
Averaging over spin directions nullifies the spin terms of the interaction potential. 
The so-called cut-off terms of the Bonn potential can be neglected, since they start 
playing an essential role only for very short distances < 0.1 fm , which would 
correspond to the baryon density % > 10 3 n # • We assume that the interaction 
between any pair of nucleons is the same, because of which in the isospin term of 
the Bonn potential we put the total isospin I\ + 1 2 — 1 describing the interaction 
between protons or neutrons. The so obtained radial part of the Bonn potential 
[113] reads 

1 Oti 



*33(r)=E-exp(-7,r) (48) 
i=i r 

with the parameters 

cti = 16.7, a 2 = 2.7, a 3 = —7.8, a 4 = —2.7, 

71 = 738 MeV, 7 2 = 769 MeV, 73 = 550 MeV, 74 = 983 MeV 

The interaction potential (48) is integrable, thus, it does not necessarily require the 
smoothing procedure [35,88,114]. For the interaction-energy density (29) we have 

$33 = / <5>3 3 (r)d7, (49) 

which yields $ 33 = 4.1 x 10" 5 MeV~ 2 = 315 MeV fm 3 . Note that $ 33 po = 

]/g 2/3 

164 MeV , hence $33Po ~ E = p ' , from where $33 ~ p . 
For the spectra of baryons we take 

£i (k) = ^k 2 + m 2 + I (p - p q ) $33, (50) 
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where p q is the quark density and i enumerates multiquark clusters: 3g, 6q, 9q 
12q , and so on. The masses of bound clusters up to z\ = 12 are taken from Table 
1, with the six-quark mass m 6 = 1944 MeV . For Zi > 15 we use the formula 
rrii ps (zi/3)m 3 for the masses of heavy multiquark clusters [22]. 

For quarks we accept the mass m q = 7 MeV and the degeneracy factor 
( q = 12 corresponding to spin 1/2, N c = 3 , and Nf = 2 . The bag constant 
B l ' A = 235 MeV . 

The results of numerical calculations [115-119] are presented in Figs. 6-13. At 
= and p = Pq , the 6q -probability is Wq = 0.18 , which agrees with the 
estimates of the 6q -admixture in nuclei [120] . The heavy-multiquark probabilities 
are always small: w$ < 10~ 3 , Wu < 10~ 5 , and W15 < 10~ 7 . At zero temperature, 
only the Bose-condensed 6q -clusters exist, the probabilities of heavier ones being 
strictly zero. Unbound quarks, at = 0, are absent below the density p n q uc ps 2p 
when they start appearing. This is why the characteristic density p q uc may be 
called the nucleation density. As we see, the value of the latter agrees with the 
corresponding estimates from Introduction. 

The stability of the quark-baryon mixture is controlled by checking the minimum 
of the free energy F — Q + J2i ^i^i whose density can be written as 

F n 

J = y = y = ~P + W, 

i 

and also by requiring the validity of the stability conditions [54] 

The probability of unbound quarks increases, with temperature or density mono- 
tonically showing that the deconfmement is a gradual crossover but not a sharp 
transition. This is in agreement with numerical simulations [121] on 16 3 x 24 
lattice, for Nf = 2 , which has demonstrated that quark-quark correlation func- 
tions at ps 1.50^ are very similar to the zero-temperature wave functions of the 
corresponding particles. 



12 Zero Baryon Density 

The case opposite to that of the previous section is when the baryon density is zero, 
n 0B = , and all particles are generated from vacuum. Then ^ — . In this 
case, we can compare our calculations with the lattice numerical simulations that 
are available only for n B = 

The spectra of gluons and quarks are again taken in the form 

B B 

uj g (k) = k + -, uo q {k) = ^k 2 + m 2 q + - (51) 

with the bag-motivated mean fields. The interaction of hadrons is considered in the 
excluded- volume approximation. The cluster volumes Vi = mjt> 2 /m 2 , according to 
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(45), are expressed through the volume v-i = 47rr|/3 of the lightest cluster with 
Zi — 2 . The bag constants for the SU (2) and SU (3) systems are to be different 
[122] with the relation B S u{2) ~ 0AB S u(3) ■ 

The presentation of results is convenient to perform in relative quantities reduced 
to these of a reference system. The role of such a system is naturally played by the 
Stefan-Boltzmann quark-gluon gas. The latter, by definition, is an ensemble of free 
massless quarks, antiquarks, and gluons. The pressure and energy density of the 
Stefan-Boltzmann quark-gluon plasma are, respectively, 

Psb = pf + pf + P f\ e SB = ef + 4 0) + ef, (52) 

where 



pf ] = ±QQ J In [l±nf\k) 



d k 



(2tt) 3 ' 

ef^QjknPik)^ (53) 

the index i — q,q,g enumerates quarks, antiquarks, and gluons with the momentum 
distribution 

nf )) (k) = {exp[P(k- f x l )]Tir 1 , (54) 

in which the upper sign is for Bosons (gluons) and the lower one, for Fermions 
(quarks and antiquarks); the chemical potentials being 

= -Vq = V, ^9 = °- ( 55 ) 

Eqs. (54) and (55) permit to write down (53) as 

,(o)_ c< r k3dk ,(o)_,jo) 



Pl ~6Wo exp|)9(ife-^)]=Fl' £l ~ 3Pl ' (56) 



An exact integration yields 

P q +Pq - 12 1 30 + ^ + 2tt 2 



(0) , „(0) _ Cq ( ^4 , ij2c> 2 |_ I 1 



7T 2 



Vf = ^ (57) 
Thus, the Stefan-Boltzmann pressure is 

^-s(^^) 94 4" !e2 ( 1+ w)- < 58 > 

This is to be compared with the QCD pressure corresponding to the QCD grand 
potential, given in Section 3, for zero coupling g — , 

P Q co = g(^"l + \n,n) S* + (l + ^) . (59) 
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Due to the relations for the degeneracy factor of gluons, ( g = 2 x (A" 2 — 1) , and 
for that of quarks, Q = 2 x Nf x N c , and antiquaries, ( q = Q , we see that (58) 
and (59) coincide with each other. Therefore, the Stefan-Boltzmann plasma is the 
asymptotic high-temperature limit of quantum chromo dynamics. 
The baryon density for the Stefan-Boltzmann plasma is 

n B = \ (Pq ~ Pi) = J J M k ) ~ U i k )\ ( 60 ) 

Either calculating (60) directly or using the derivative % = dp/d/iB , with /i B = 
3/x we have 

nB = C «lfe v" 2 + 7r2@2 ) • (61) 
From here, one gets the equation 

^3 + 7r 2 Q 2 /i _ 18^ nB = q (62) 

defining n = [1(113) ■ At zero baryon density n# = , as is clear from (62), one 
has fi = . 

When the chemical potential is zero, the density of quarks becomes 

= C*4^fC(3) (0 = 0), (63) 

where £(3) = 1.20206 . For gluons, the chemical potential is always zero, so their 
density is 

P9 = C^C(3). (64) 
Finally, for the specific heat of the Stefan-Boltzmann plasma we find 

de SB 2tt 2 / 7 \ 3 1 2m2 ^ 2 - 7T 2 2 
^ = - 50- = 15 { C ° + 9 + 2 C ^ 6 VT^ ' (65) 

The results of numerical calculations [123-125] for the mixed system, consisting 
of the quark-gluon plasma with the spectrum (51) and of hadrons in the excluded- 
volume approximation, will be presented below for three different situations. 

12.1 SU(2) Quarkless System 

The system consists of unbound gluons and of glueballs that are bound gluon clus- 
ters. The experimental status of glueballs is yet uncertain, though there are sugges- 
tions [126,127] to interpret a narrow resonance appearing in proton-proton collisions 
as a scalar glueball. Pure gluodynamics is often studied because it is easier, than 
the full chromodynamics, for Monte Carlo lattice simulations. 

Glueball masses have been computed in the lattice gauge theory for both SU (2) 
[43,46] as well as for SU(3) [128-130] cases. The lattice results are close to the 
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bag-model calculations [131-133]. Here we accept the glueball masses found in the 
bag-model approach [132,133]. The corresponding glueball characteristics are given 
in Table 2. The radius r 2 of the lightest glueball with m 2 = 960 MeV is a fitting 
parameter which is taken as r 2 = 1.2 fm . The constant B in (51) is chosen to be 
B = (165 MeV) A . The gluon degeneracy factor is ( g = 6 for the SU{2) case. 

The results of our calculations are displayed in Figs. 14-18, where the glueball 
probability w G and the gluon probability w g are defined by 

qlueballs 

_ Pi _ Pg n 

w G = > z~, w g = — = 1 - w G . 

i P P 

The relative energy density is compared with the lattice data [134,135]. The refer- 
ence Stefan-Boltzmann plasma here corresponds also to the quarkless case, Nf — . 
Deconfinement occurs at = 215 MeV as a second-order transition, which is in 
agreement with the lattice simulations. 

12.2 SU(3) Quarkless System 

The glueball parameters are taken from Table 2. The constant B in spectra (51) 
is B = (235 MeV) 4 . The gluon degeneracy factor for the SU(3) case is Q — 16 . 

Varying the radius r 2 of the lightest glueball, we have three possibilities: (i) 
r 2 < r c , where r c = 0.8/m ; then deconfinement is a gradual crossover, (ii) 
r 2 = r c ; in this case deconfinement is a 2 -order transition, (iii) r 2 > r c ; then 
1 -order transition occurs. These possibilities are illustrated in Figs. 19-23, where 
the entropy density at Hi — is s = /3(£+p) and the reference Stefan-Boltzmann 
entropy density is ssb = P(£sb + Psb) = 4=Ppsb ■ The relative entropy density is 
compared with the lattice numerical simulations [136]. The latter agrees with our 
results if r 2 > r c , so that deconfinement becomes a 1 -order transition at about 
B d = 230 MeV . 

Emphasize the importance of taking into consideration glueball interactions: 
When these are absent, that is r 2 = , the behaviour of the system is unphys- 
ical. 

12.3 SU(3) System with Quarks 

The constituents of the system are taken as follows. Consider quarks of two flavours, 
q = {u, d} , and the corresponding antiquarks q = {u, d} with the masses m q = 
rriq = 7 MeV . The degeneracy factor for each pair of up and down quarks is 
( q = 2 x Nf x N c = 12 , and the same for antiquarks, (g = 12 . Gluons have 
the degeneracy factor ( g = 2 x (JV^ — 1) = 16 . From the long list of the known 
hadrons, we include only those with the lightest masses, which mainly contribute to 
thermodynamics. These are unflavoured mesons (Table 3), strange mesons (Table 
4) and light baryons (Table 5). 

For the radius of the lightest hadron, that is of a pion, we take r 2 = r n = 
0.56 fm . The radii of all other clusters are expresses through r n using (45). For 
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the mean-field parameter in (51), we accept B 1 ^ = 210 MeV . The results of 
calculations are shown in Figs. 24-26, where the hadron cluster probability is defined 

as 

clusters 

\ ^ Pi 

W c = > Zi — . 

i P 

Deconfinement is found to be rather a continuous crossover-like transition at 0^ = 
166 MeV , which is close to lattice data [137]. 



13 Thermodynamic Restriction Rule 

Invoking this or that approximation, one gets an effective thermodynamic potential, 
for instance, an effective grand potential Q = Q(0, V, p, tp) , involving some auxil- 
iary functions tp = {<Pj} depending on thermodynamics parameters, temperature 
O , volume V , and a set p = {p{\ of chemical potentials. Thus, effective spectra 
in (51) contain the mean field tp = B / p . In the excluded-volume approximation, 
the free volume of the system is factored with the quantity £ = 1 — J2i Pi v i > as is 
seen from (42). Both p and pi are functions of 0, V, p . In the cut-off model of 
Section 5, the effective cut-off momentum k is a function of . 

Effective thermodynamic potentials are to be handled with great caution. Really, 
if one calculates, e.g., the pressure 

p = -dv = -v = v lnTie (66) 

in two different ways, as the derivative (— dVt/dV) or as the ratio (— Q/V) , then 
the answers can be different when Q includes auxiliary functions depending on 
V . This would mean that the relation (66) breaks. The same concerns the energy 
density, entropy density, and the cluster densities, respectively, 

dv ' 



dp 



f3^£ + P-J2 PiP^j > 



Pi * 

which may be defined in two ways, as first derivatives of pressure or as the corre- 
sponding statistical averages. Definition (67) can also become broken when auxiliary 
functions depend on thermodynamic variables. This kind of inconsistency occurs as 
well for the second derivatives, such as the specific heat 

c " = § = v ( {g2} ~ (68) 
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or the isothermic compressibility 

4 «*>-<*>')■ (69) 

These inconsistencies in defining thermodynamic characteristics in two ways, 
thermodynamic and statistical, of course, are not pleasant. Moreover, the difference 
between these two ways is not only quantitative, but can also become drastic, espe- 
cially for systems with phase transitions. It is even possible to give examples when 
the definition through the derivatives yields unphysical divergencies in the energy 
and entropy densities at the phase transition point. This, for instance, happens, as 
is easy to check, for a pure gluon model in the effective spectrum approximation. 

The simplest procedure for avoiding the described troubles can be formulated 
as follows. Let x be any of the thermodynamic variables 0, V or /j, . If 
Q is an effective grand potential including auxiliary functions depending on these 
thermodynamic variables, then 

dtt _ fdfl\ Oil dip 
dx \ dx J dp dx 

It is just the second term here which causes all unpleasant problems. So, the decision 
is evident: the derivatives dfl/dx are to be understood in the restricted sense as 

dVl ^ fdtt\ 
dx \dx J 

The consent (70) may be called the thermodynamic restriction rule. We always 
employ this rule dealing with effective thermodynamic potentials. If the derivatives 
in (66)- (69) are understood in the sense of (70), then both ways of calculating 
thermodynamic characteristics yield the same answers. 

Although with the restriction rule (70) we avoid the appearance of spurious 
terms so that all relations (66)-(69) become self-consistent, another problem can 
arise when dealing with effective thermodynamic potentials. This is the occurrence 
of instability regions around a transition point, where either the specific heat Cy 
or the isothermic compressibility kt are negative. Below we illustrate this for 
the SU(3) quarkless system of subsection 12.2 with the mean-field parameter 
B = (225 MeV) 4 . The results are shown in Figs.27-34. The unstable solutions 
appearing in the vicinity of transition points are related to the loss of convexity 
of pressure. To restore the convexity, we may resort to the Maxwell construction 
smoothing the corresponding thermodynamic potential [138]. Such a smoothing 
is shown in Figs. 31 and 32. The behaviour of the resulting pressure and energy 
density is in a reasonable agreement with lattice simulations [50]. Nevertheless, 
a slight dissatisfaction rests with the fact that instability regions occur not only 
around a 1 -order transition, where this would be more or less natural, but also in 
the vicinity of a continuous transition. 

Instead of relying on the restriction rule (70), it would seem rational to define 
an effective thermodynamic potential, from the beginning, in such a way that all 
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(70) 



necessary thermodynamic relations would be automatically valid. This goal can 
be achieved [139] by redefining thermodynamic characteristics with the help of the 
shifts of the chemical potentials, //j — > fii — Ui , pressure, p — > p+p' , energy density, 
e — > e + e' , and entropy density, s — > s + s' requiring that the shifting functions 
Ui,p',e' , and s' guarantee the validity of (66) and (67). The latter then are 
named the self-consistency conditions [139]. In this case, (66) and (67) is a system 
of nonlinear differential equations, for the functions Ui,p',e' and s' , in partial 
derivatives with respect to the variables 6, V and ^ . Such a system has no unique 
solution, especially when boundary conditions are not known. To extract a solution 
from the self-consistency equations needs several additional heuristic assumptions 
and fitting parameters. Some simplification comes from the guideline prescribed by 
mean-field approximations [140-142]. 

14 Principle of Statistical Correctness 

In this section we present a new principle allowing a correct construction of effective 
thermodynamic potentials. This principle, as compared to the self-consistency con- 
ditions, is: (i) more general, yielding these conditions but not conversely; (ii) much 
simpler to deal with; (iii) unambiguous, providing a unique solution. 

Let an effective thermodynamic potential fi e // = ^effif) include a set ip — 
{{pi(x)} of auxiliary functions depending on arbitrary variables x . The latter, in 
particular, may incorporate space and thermodynamic variables. First of all, it is 
necessary to understand that not any effective potential can have sense, however 
reasonable it may look. Each thermodynamic potential, to be accepted as such, 
must have the properties formulated below. 

Property 1. Statistical Representability: 

An effective thermodynamic potential Q e ff represents an equilibrium statistical 
system if and only if it has the Gibbs form 

n eff (< P ) = n[H eff (<p)], (7i) 

where 

tt[H] = -GhiTr e~ m , (72) 

depending on auxiliary functions only through an effective Hamiltonian H e ff = 
H e ff(ip) . Such a thermodynamic potential is called statistically representable. 

In this way, if one invents an effective thermodynamic potential, even pronounc- 
ing seemingly plausible words, this does not mean that the invented potential de- 
scribes some statistical system. If the potential is not statistically representable it 
represents no equilibrium statistical system. For example, a thermodynamic po- 
tential in the excluded-volume approximation is not statistically representable. Al- 
though the latter approximation may occasionally give a reasonable description, 
but in general it is not trustworthy. The excluded-volume approximation may be 
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used, because of its simplicity, as a first attempt of understanding the qualitative 
behaviour of a system, but it should be always followed by a more reliable approxi- 
mation. 



Property 2. Thermodynamic Equivalence: 

A statistical system described by a given Hamiltonian H g i V is thermodynami- 
cally equivalent to a system modeled by an effective Hamiltonian H e ff if and only 
if their thermodynamic potentials are statistically representable, 

Q giv = Q[H giv ], Q e ff = Q[H e ff], (73) 

and are equal to each other, 

n[H giv ] = n[H eff \. (74) 

The corresponding Hamiltonians are called thermodynamically equivalent. 

For the case of infinite matter, such as nuclear matter, the equality (74) can be 
softened by requiring the validity of the asymptotic, in the thermodynamic limit, 
equality 

lim 1 (Q[H giv \ - Q[H eff \) = 0. 



Property 3. Statistical Equilibrium: 

The necessary condition for an equilibrium statistical system modelled by an 
effective Hamiltonian H e ff(<p) to be thermodynamically equivalent to a given 
statistical system with H giv is the equilibrium condition 

where the variation over ud implies a set of variations with respect to each udi and 

~ TrAexp(-pH eff ) 
[ ' Trexp(-(3H eff ) ' 

The proof of (75) is straightforward basing on the statistical representability 
(71), thermodynamic equivalence (74) and the fact that Q g i V does not depend on 
up . 

Now we can formulate the central notion: 
Principle of Statistical Correctness: 

An effective thermodynamic potential is statistically correct if it is statistically 
representable with an effective Hamiltonian satisfying the condition (75) of statistical 
equilibrium. 

As is evident, the self-consistency conditions for the first-order derivatives (66) 
and (67) immediately follow from (75). Moreover, the self-consistency conditions 
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for the second-order derivatives (68) and (69) also follow from (75) as well as such 
conditions for the derivatives of arbitrary order. While if one finds the shifting 
functions from the first-order self-consistency conditions (66) and (67), the second- 
order conditions (68) and (69) are not necessarily fulfilled. 

We shall also say that an effective Hamiltonian is statistically correct if it satisfies 
(75). The same can be said about an approximation leading to a statistically correct 
Hamiltonian. For instance, the correlated mean-field approximation of Section 9, 
involving conditions (26) or (27), or (31), is statistically correct. 



15 Clustering Quark-Hadron Matter 

To obtain a statistically correct description of the quark-gluon plasma clustering 
into hadron states, let us use the correlated mean-field approximation [35] leading 
to the Hamiltonian 

H = J2H t + CV, 

i 



LO, 



i(k) = y/tf + mZ + Ui-fMi (76) 

discussed in Section 9. 

Consider the case of the conserved baryon number Nb = J2i with Nf = 
BiNi , where Bi is the baryon number of an % -cluster. For an equilibrium system, 
the relation 

Hi = BifiB (77) 

holds between the chemical potential pi and the baryon potential ps ■ The latter 
may be defined as a function psins) of the baryon density 

n B = ^- = Y, B iPi- (78) 

i 

The index % enumerates the constituents. The total set {i} of these indices 
consists of two different groups, {i} = {i} p i U Wc/ • The first group {i} P i 
corresponds to the plasma constituents, quarks, antiquarks, and gluons, which are 
elementary particles, thus, having the compositeness number Zi — 1 . The second 
group {i} c i enumerates hadron clusters that are bound states with compositeness 
numbers > 2 . Respectively, the total density of matter 

p = J2 z iPi = P P i + Pch (79) 

i 

consists of two terms 

Ppi = S P*> Pci=Y, Z iPi 

{i} P l {i}d 

corresponding to the plasma density p p i and the cluster density p c \ . 
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Define the plasma mean fields Ui , when % e {i} p i , as 

Ui = U(p) =pj V{r)s{r)d r, (80) 

where V(r) is a confining potential and s(r) , screening function. Before substitut- 
ing into (80) a concrete confining potential, let us emphasize the general properties 
which the plasma mean field U (p) must satisfy to. These properties are 

l7(p)_oo (p_0), 

U(p) ^ (p -> oo). (81) 

The upper line in (81) means that quarks and gluons cannot exist as unbound 
particles at low density, that is, the colour confinement must occur as p — > . Or one 
may say that quarks and gluons cannot exist as free particles outside dense nuclear 
matter. The lower line in (81) reflects the phenomenon of asymptotic freedom. 

There are different types of confining potentials, linear, quadratic, logarithmic, 
and with noninteger powers. For example, the interaction between a heavy quark 
and its antiquark is usually taken in the form of the Cornell potential [143,144] with 
the linear confining term. This form of the potential is confirmed by QCD calcula- 
tions [145] and by lattice simulations [146]. The quadratic confining potential is also 
quite popular [77,80,91]. In addition, the intensity of mutual interactions between 
three plasma constituents, quarks, antiquarks, and gluons, is different [145,147]. The 
confining potential V(r) in (80) is assumed to be an averaged potential of the form 

V(r) = Ar v (0 < v < 2). (82) 

The screening function s(r) = c(r/a) is usually [91] scaled with the mean 
interparticle distance a = p -1 / 3 . Therefore, the plasma mean field (80) with the 
confining potential (82) can be written as 

U(p) = J^V 173 , (83) 
where _ 

POO 

J 1+u = An A / c(x)x 2+u dx. 
Jo 

We can calculate the constant J if A and c(x) are known. Alternatively, we 
may treat J as a free parameter. The value of J can be estimated as follows. 
Accept that at the normal quark density p the plasma mean field (83) becomes 

U(p ) = 3E = 3pl /3 , (84) 

where the factor 3 stands for the three plasma constituents. Then from (83) and 
(84) we obtain 

J = ^ 1+ ^pl'\ (85) 

Thus, for the linear confinement, v — 1 , we get J = 272 MeV , while for the 
quadratic confinement, v = 2 , we have J = 226 MeV . 
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For the mean field of an i -cluster we may write 

Ui = £ ®aPj + *i P(p) - U(p d )} , (86) 

{j}cl 

where the first term describes the interaction of the given cluster with other clusters, 
and the second term corresponds to the interaction of this cluster with the quark- 
gluon plasma. The interaction potentials between clusters can be scaled according 
to (20), which permits to express the interaction integrals (29) through one scaling 
integral as 

= ZiZj®. (87) 
Taking into account (83) and (87) reduces (86) to 

Ui = z^p d + ZiJ 1+v [p- vjz - Pc 7 /3 ) . (88) 

In the effective Hamiltonian (76) with the mean fields (83) and (88) the role of 
auxiliary functions is played by the densities p and p c i . So, for the equilibrium 
conditions (75) we have 

/9H\ 

(W=°- (89) 

From (89) we have two variational equations of the type (31), whose solution, up to 
a constant, is easy to find: 

c = (V-^ 3 - P ir /3 ) - \*& (90) 

In this way, the effective Hamiltonian is completely defined and we can pass to 
particular calculations. 

15.1 Pure Gluon System 

Imagine an extreme situation when only gluons can exist. This case may be obtained 
from the general model by putting all degeneracy factors zero except that of gluons, 
Q 7^ . Then p = p g and p c i — . Employing the Boltzmann approximation, we 
find [148] that there occurs a first order transition vacuum-gluon plasma at 




<dd = J 



3-i/ 



exp 




V(i+") 



Below there is exactly vacuum, empty space, with zero energy density 5 = 0. 
Gluons appear at 0^ with a jump. The relative latent heat at 0^ is 

Ae d 1 + u A 3 s 



exp 1 

The degeneracy of gluons for the SU(3) case is (=16. For the linear comfinement 
with v — 1 , we get 0^ = 248 MeV and Aea/esB = 0.27 . For the harmonic 
confinement, when v — 2 , the vacuum-gluon transition happens at a higher 
temperature 0^ = 285 MeV and the latent heat is larger, Aed/ssB = 0.91 . 
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15.2 SU(2) Gluon-Glueball System 

The gluon degeneracy factor for the SU (2) case is ( = 6 . As the scaling integral we 
take $ = $22/4 . So that the interactions between glueballs are found from (87). We 
consider the glueballs listed in Table 2. There are two fitting parameters for which 
we accept J = 175 MeV and $ 22 = 38.42 GeV fm 3 , so that $ = 9.61 GeV fm 3 . 
The results of our calculations [148-151] are shown in Figs. 35-37 for the quadratic 
confinement, v m 2 . Actually, the results do not change much in the interval 
1.5 < v < 2 . We have mainly used v = 1.86 . Deconfinement is a second-order 
transition at Q d = 210 MeV . As is seen, the agreement with the lattice simulations 
[135] is beautiful. 

15.3 SU(3) Gluon-Glueball Mixture 

For the gluon degeneracy factor we have ( g — 16 . The power v of the confin- 
ing interaction is, as in the previous subsection, close to quadratic. But the fitting 
parameters are J = 225 MeV and <£> = 3.84 GeV fm 3 . The glueball characteris- 
tics are again taken from the Table 2. Calculations show [148,149] that a first-order 
transition occurs at 0^ = 225 MeV with the relative latent heat Ae d /e 5B = 0.23 . 
The agreement with the Monte Carlo lattice simulations [50,51,152] is also very good. 

15.4 Mixed Quark-Gluon-Meson System 

Consider the case of zero baryon density, % = . Take the quarks of two flavours, 
u and d , and the related antiquarks, u and d . Assume, for simplicity, that 
all their masses are equal, m u = rrid = 7 MeV . The degeneracy factor for each 
kind of quarks is ( u = (d = 2A" C = 6 . This factor for gluons is ( g = 16 . Hadrons 
are represented by mesons from Tables 3 and 4. Following Section 11, we accept 
$33 = 315 MeV fm 3 , so that $ = $ 33 /9 = 35 MeV fm 3 . The plasma interaction 
parameter J = 225 MeV is the same as in the previous subsection, as well as 
v ?s 2 corresponding to quadratic confinement. So, here we do not add any new 
fitting parameters. 

Our calculations [148,149] displayed in Fig. 40 prove that there is no sharp phase 
transition but there is a gradual crossover. The deconfinement transition can be 
attributed to the temperature <d d = 150 MeV where the relative specific heat 
Cy /Csb has a maximum. The latter is finite and look rather as a Schottky anomaly 
[153] than as a narrow divergent peak typical of a second-order phase transition. 
The agreement of our results with the lattice-simulation data [137] is again quite 
good. The lattice results [137] indicate that the deconfinement transition is really 
continuous. From the point of view of QCD this can be understood as follows. The 
role of the quark term in the QCD Lagrangian is similar to that of an external 
magnetic field applied to a spin system. In the presence of a magnetic field, the 
ferromagnet-paramagnet transition in simple spin systems becomes a continuous 
crossover. 
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15.5 Finite Baryon Density 



Here we extend the consideration to nonzero baryon density. The parameters J 
and $ are the same as in the previous subsection. Again, we study the two- 
flavour case with the same characteristics. We take mesons from Table 3, protons 
and neutrons from Table 5, and multiquarks from Table 1. For a six-quark cluster 
we accept m 6 = 1944 MeV . The Bose-Einstein condensate of six-quarks occurs 
when ^(O) = . Then, the baryon potential is 

\ib = \me + 3$p ci + 3 J l+V (p^ 3 + p c 7 /3 ) , 
and the density of six-quarks consists of two terms: 

Pe = J n 6 (k)d k +p° 6 . 

We have analysed the behaviour of several probabilities [154-157] as functions 
of temperature 6 and relative baryon density n B /n B ■ This is demonstrated in 
Fig. 41 for the plasma probability 

Wpl = ~(Pg + Pu + Pu + Pd + Pd) , 

in Fig. 42 for the pion probability 

2 / 

W n = ~ (Pk+ + Pit- + P*°) , 
P 

in Fig. 43 for a summarized, excluding pions, probability of other mesons 

2 / 

W V pu = ~ [P V + Pp+ + Pp- + Pp" + Pu>) = W mes , 

in Fig. 44 for the nucleon probability 

3 / 

W 3 = - (Pp + Pp + Pn+ Pn) , 
P 

in Fig. 45 for the six-quark probability 

6 < 

w 6 = - (p 6 + pe) , 
P 

and in Fig. 46 for the probability of condensed six-quark clusters 

- 6 ° = - p pl ■ 

In addition, we present here some other thermodynamic characteristics permit- 
ting to better understand the features of the deconfinement transition. The ratio of 
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pressure over energy density, which has the meaning of the effective sound velocity 
squared, 

c 2 - V - 

is given in Fig. 47. At n# = the temperature dependence of c 2 e ^ agrees with 
that reconstructed from the lattice data [158]. The reduced specific heat 



are depicted in Figs. 48 and 49, respectively. The transition line can be ascribed 
to the maximum of the inverse compressibility coefficient k^. 1 which can be called 
[111] the compression modulus (see Fig.50). 

We shall not discuss in detail the peculiarities of the calculated thermodynamic 
functions. This is because, as we think, the presented figures already give a good 
visual demonstration, and also in order not to make this review too long. Let us 
only emphasize that the deconfinement transition is a continuous crossover becoming 
smoother and smoother with increasing baryon density. The deconfinement at a 
fixed low temperature and rising baryon density is due to the disintegration of 
hadrons into unbound quarks. When both temperature and baryon density increase, 
the deconfinement is a result of the hadron disintegration as well as of the generation 
from vacuum of quarks and gluons. 

Concluding we may state that taking into account the coexistence of hadrons and 
of the quark-gluon plasma is vitally important for constructing a unified approach 
being in agreement with the lattice-simulation data. The gradual character of the 
deconfinement transition, occurring through a mixed hadron-plasma state, rules 
out those predictions that have been based on a sharp first-order phase transition. 
This concerns the interpretation of signals of the quark-gluon plasma at heavy ion 
collisions [8-14,159] and the hadronization scenario related to the evolution of early 
universe after the Big Bang [4]. The quantitative predictions of our approach can 
be improved in several ways, for instance, by including more kinds of particles or 
by invoking more elaborate interaction potentials [160]. However, we do hope that 
qualitatively the picture will remain the same. 
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Figure Captions 
Fig.l. 

Illustration of a cluster fusion. 
Fig.2. 

6q -probability vs. the mass of a 6q -cluster at 6 = and p = p . 
Fig.3. 

Cluster probabilities as functions of the relative density at = 200 MeV . 
Fig.4. 

Nucleon probability as a function of the relative density at different temperatures. 
Fig.5. 

6q -probability vs. relative density at = 0. 
Fig.6. 

Nucleon, 6q -cluster, and unbound quark probabilities as functions of the relative 
density at = 0. 

Fig. 7. 

Probabilities of unbound quarks and of bound clusters vs. relative density at = 
30 MeV . 

Fig.8. 

Cluster probabilities vs. relative density at = 50 MeV . 
Fig.9. 

Cluster probabilities at = 100 MeV . 
Fig.10. 

Cluster probabilities at = 200 MeV . Dashed lines show the points of first order 
phase transitions. Between these points the matter is a stratified gas-liquid mixture. 

Fig.ll. 

Nucleon, 6q -cluster, and quark probabilities as functions of temperature in MeV 
at the normal density p = po . 

Fig.12. 

Cluster probabilities vs. temperature at the fixed density p = 5p • 
Fig.13. 

Phase portrait for the baryon rich quark-hadron matter. Along the dashed line the 
compressibility is divergent. 
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Fig.14. 

Glueball probability for the SU (2) quarkless system as a function of temperature. 
Fig.15. 

Comparison of glueball and gluon probabilities for the SU (2) quarkless system. 
Fig.16. 

Relative energy density for the SU(2) quarkless system (solid line) as compared 
with the lattice Monte Carlo data (Engels et al., 1981). 

Fig.17. 

Relative energy density for the SU (2) quarkless system (solid line) compared with 
the lattice numerical simulations (Engels et al., 1989). 

Fig.18. 

Reduced specific heat for the SU (2) gluon-glueball mixture. At the deconfinement 
temperature, specific heat diverges. 

Fig.19. 

Glueball probability for the 577(3) quarkless system at several values of the lightest 
glueball radius: (1) r 2 = ; (2) r 2 = 0.5 fm ; (3) r 2 = 0.7 fm ; (4) r 2 = 0.8 fm . 

Fig.20. 

Relative energy density for the S77(3) quarkless system at: (1) r 2 = ; (2) 
r 2 = 0.5 fm ; (3) r 2 = 0.7 fm ; (4) r 2 = 0.8 fm ; (5) r 2 = 1 fm . 

Fig.21. 

Relative entropy density for the SU(3) quarkless system at r 2 = 0.82 fm 
compared with the lattice numerical data (Brown et al., 1988). 

Fig.22. 

Reduced specific heat for the S77(3) gluon-glueball mixture at r 2 = 0.8 fm . 
Fig.23. 

Relative energy density for the SU (3) quarkless system at 77 1 / 4 = 210 MeV and 
r 2 = 0.6 fm . 

Fig.24. 

Hadron cluster probability for the mixture of quark-gluon plasma and hadrons at 
zero baryon density. 

Fig.25. 

Relative energy density for the mixture compared with the lattice numerical calcu- 
lations (Qelik et al, 1985). 
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Fig.26. 

Reduced specific heat for the mixture of quark-gluon plasma and hadrons at zero 
baryon density. 

Fig.27. 

Total density as a function of temperature for the SU (3) gluon-glueball mixture at 
different values of the lightest glueball radius: (1) r 2 = 0.6 fm ; (2) r 2 = 0.8 fm ; 
(3) r 2 = 1 fm . 

Fig.28. 

Gluon probability vs. temperature for the values of r 2 as in Fig.27. 
Fig.29. 

Pressure of the gluon-glueball mixture vs. temperature for the same values of r 2 
as in Fig.27. 

Fig.30. 

Relative energy of the mixture for the values of r 2 as in Fig.27. 
Fig.31. 

Smoothing of pressure in the unstable crossover region at r 2 = 0.7 fm . 
Fig.32. 

Smoothing of pressure in the region of first-order phase transition at r 2 = 1 fm . 
Fig.33. 

Relative smoothed pressure for the gluon-glueball system at r 2 = 0.9 fm (solid 
line) as compared with the lattice simulation data (Engels, 1991). 

Fig.34. 

Relative smoothed energy at r 2 = 0.9 fm (solid line) compared with the lattice 
data (Engels, 1991). 

Fig.35. 

Gluon probability for the corrected SU (2) quarkless model. 
Fig.36. 

Relative energy and pressure for the corrected SU(2) gluon-glueball model. Circles 
and squares are lattice simulation data (Engels, 1989). 

Fig.37. 

Reduced specific heat for the corrected SU (2) quarkless model. 
Fig.38. 

Gluon probability for the corrected SU (3) gluon-glueball model. 
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Fig.39. 

Relative energy and pressure for the corrected 577(3) quarkless model compared 
with the lattice Monte Carlo calculations (Brown et al., 1988; Karsch, 1989; Engels, 
1991; Petersson, 1991). 

Fig.40. 

Relative energy and pressure for the mixed quark-gluon-meson system. 
Fig.41. 

Quark-gluon plasma probability on the temperature-baryon density plane. 
Fig.42. 

The probability of ir -mesons. 
Fig.43. 

The total probability of 77 -, p -, and uj -mesons. 

Fig.44. 

Nucleon probability. 

Fig.45. 

The probability of six-quark clusters. 
Fig.46. 

The probability of condensed six-quark clusters. 
Fig.47. 

Pressure-to-energy ratio on the temperature-baryon density plane. 

Fig.48. 

Reduced specific heat. 

Fig.49. 

Dimensionless compressibility coefficient. 
Fig.50. 

Dimensionless compression modulus as a function of temperature and baryon den- 
sity. 



51 



Table Captions 
Table 1. Multiquark parameters 
Table 2. Glueball parameters 
Table 3. Unflavoured meson parameters 
Table 4. Strange meson parameters 
Table 5. Baryon parameters 
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Table 1 



mass 


compositeness 


degeneracy 


mi (MeV) 


number z\ 


factor Q 


939 


3 


4 


1944 


6 


9 


2163 


6 


3 


3521 


9 


4 


4932 


12 


1 



Table 2 



mass 


compositeness 


degeneracy 


m { (MeV) 


number Zi 


factor Q 


960 


2 


6 


1290 


2 


6 


1590 


2 


6 


1460 


3 


11 


1800 


3 


39 



Table 3 



mesons 


mass 


compositeness 


degeneracy 




mi (MeV) 


number Zi 


factor Q 


7T + 


140 


2 


1 


7T~ 


140 


2 


1 


7T° 


135 


2 


1 


V 


548 


2 


1 


p + 


770 


2 


3 


p 


770 


2 


3 


p° 


770 


2 


3 


UJ 


782 


2 


3 
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Table 4 



mesons 


mass 


compositeness 


degeneracy 




mi (MeV) 


number Zi 


factor Ci 


K + 


494 


2 


1 


K~ 


494 


2 


1 


K° 


498 


2 


2 


K° 


498 


2 


2 



Table 5 



baryons 


mass 


compositeness 


degeneracy 




mi (MeV) 


number z-i 


factor Q 


N 


939 


3 


4 


N 


939 


3 


4 


A 


1232 


3 


16 


A 


1232 


3 


16 
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